#
# MeasureMisperceptionAnalysis.R
#
# Replication file for Carlson, Taylor N. and Seth J. Hill, 
# ``Experimental Measurement of Misperception in Political Beliefs.''
# Journal of Experimental Political Science.
#

rm(list=ls())
if (.Platform$OS == "windows") {
  # Set working directory in location of script.
  .doit <- function() { # only works with R.exe; trap errors if using Rscript.exe
  f_f <- lapply(sys.frames(),function(x) x$ofile);f_f <- Filter(Negate(is.null),f_f) ; PTH <- dirname(f_f[[length(f_f)]]); setwd(PTH) ; rm(PTH,f_f)
  }
  try(.doit(),silent=T)
}
library(bit64)
library(data.table)
options(stringsAsFactors=F)
library(RColorBrewer)
palette(brewer.pal(n=9,name="Set1")[c(1:5,7:9)]) # reset default colors
options(digits=4)

# File location for tables and figures.
.fig.path = c(".")
.tab.path = c(".")
postpend = ""

#
# Call in data and create functions.
# 
DT = fread("Data.csv")

# Recodes.
DT[,subject.income.cat := factor(subject.income.cat,levels=c("Less than $25K","$25-$40K","$40-$60K","$60-$90K","$90-$150K","$150-$250K","Above $250K"))]
DT[,income.cat := factor(income.cat,levels=c("Less than $25K","$25-$40K","$40-$60K","$60-$90K","$90-$150K","$150-$250K","Above $250K"))]
DT[,vote.16 := factor(ifelse(votechoice16 == 1,"Clinton",ifelse(votechoice16 == 2,"Trump","Other")),levels=c("Clinton","Trump","Other"))]
DT[,ub.trump.2p := ub] # previously calculated from ANES data
DT[,lb.trump.2p := lb]  
DT[,prob.in.ci := 1*(posterior.trump >= 100*lb.trump.2p & posterior.trump <= 100*ub.trump.2p)]
DT[,dist.from.share := posterior.trump - 100*trump.2p]
DT[,abs.dist.from.share := abs(posterior.trump - 100*trump.2p)]
# Distance given estimate outside of ci.
DT[,dist.from.share.notin.ci := ifelse(prob.in.ci,NA,dist.from.share)]
DT[,dist.toward.true := (posterior - prior)]


#
########################
# Functions.
########################
#
.mkCoefs <- function(lm.out,regex="full.treat") {
  # Take lm object and turn to a data.table of coefs and  ses.
  # Arg regex is regex pattern to remove from names.
  coefs = coef(lm.out)
  coefs=data.table(yname=names(coefs),x=coefs,se=sqrt(diag(vcov(lm.out))))
  coefs=coefs[regexpr("prior",yname) == -1,] # not prior
  coefs[,yname := gsub(regex,"",yname)]
  coefs[,ub := x + 1.96*se]
  coefs[,lb := x - 1.96*se]
  setkey(coefs,x)
  return(coefs)
}
.coefPlot <- function(coefs,fname=NULL,ylab="",xlab="",ynames="",cex.axis.y=.7) {
  # Make a coefficient plot and write out model names.
  # Arguments.
  #  coefs -- data.table of coefficient estimates.
  #  fname -- file name for saving figure and notes text; when NULL, plots to screen
  #  ylab,xlab -- ylab,xlab for figure
  #  ynames -- column in coefs with y-axis names.
  #  cex.axis.y -- y-axis labels cex.
  
  if (!is.null(fname)) { pdf(sprintf("../Figures/%s.pdf",fname),width=10,height=6) }
  par.old <- par(cex.lab=1.2)
  coefs[,plot(x=range(c(ub,lb)),y=c(1,.N),type='n',ann=F,axes=F)]
  grid(lwd=2);axis(1)
  abline(v=0,lty=2,lwd=2,col='gray')
  title(xlab=xlab, ylab=ylab,line=2)
  axis(2,las=2,labels=coefs[[ynames]],at=coefs[,seq_len(.N)],cex.axis=cex.axis.y)
  # CI.
  coefs[,segments(x0=lb,x1=ub,y0=seq_len(.N))]
  # Point estimate.
  coefs[,points(x=x,y=seq_len(.N),pch=19,cex=1.5)]
  par(par.old)
  if (!is.null(fname)) { dev.off() }
}
.CI95 = function(p,n,upper=T) {
  # Return upper or lower bound of CI.
  if(upper) { p + 1.96*sqrt(p*(1-p)/n) }
  else { p - 1.96*sqrt(p*(1-p)/n) }
}
.errorBars = function(z,bar.mat.ub,bar.mat.lb,lwd=2) {
  # Add error bars to a bar plot.
  # Arguments.
  #  z -- result of a barplot() call.
  #  bar.mat.ub, bar.mat.lb -- upper and lower bounds
  #   for bars.
  segments(x0=z, y0=bar.mat.ub, y1=bar.mat.lb, lwd=lwd)
  arrows(x0=z, y0=bar.mat.ub, y1=bar.mat.lb, lwd=lwd, angle=90, code=3, length=0.05)
}
.biasPlot <- function(DT,yvar="dist.from.share",xvar="prior.cat",xlab="Prior",ylab="Posterior belief minus ANES Trump share",add.lines=T,tick=T) {
  # Create the bias in perceptions plot.
  # Aggregate to mean and sd.
  bias = DT[, list("dist.mean"=mean(get(yvar),na.rm=T), "sd"=sd(get(yvar),na.rm=T),.N), keyby=xvar]
  # Append separately by characteristic just seen.
  bias = rbind(bias, DT[, list("dist.mean"=mean(get(yvar),na.rm=T), "sd"=sd(get(yvar),na.rm=T),.N), keyby=c(xvar,"input")], fill=T)
  bias[is.na(input), input := "Overall"]
  # Standard errors and 95% CI.
  bias[,se := sd/sqrt(N)]
  bias[,lb := dist.mean - 1.96*se]
  bias[,ub := dist.mean + 1.96*se]
  
  # Plot.
  bias[input == "Overall",plot(x=as.numeric(get(xvar)),y=dist.mean,pch=19,cex=1.5,type='p',ann=F,axes=F,ylim=bias[,range(c(dist.mean,0),na.rm=T)])]
  grid(lwd=2)
  abline(h=0,lty=2,lwd=3,col='gray')
  bias[,axis(1,at=seq_len(length(unique(get(xvar)))),labels=levels(get(xvar)),tick=tick)]
  axis(2,las=2)
  title(xlab=xlab,ylab=ylab)
  # Line plots for each other input.
  if (add.lines) {
    bias[input != "Overall",lines(x=as.numeric(get(xvar)),y=dist.mean,col='gray'),by="input"]
  }
  # Re-do plot for overall.
  bias[input == "Overall",points(x=as.numeric(get(xvar)),y=dist.mean,pch=19,cex=1.5)]
  if (add.lines) {
    bias[input == "Overall",lines(x=as.numeric(get(xvar)),y=dist.mean,lwd=3)]
  }
  # Error bars.
  bias[input == "Overall",.errorBars(as.numeric(get(xvar)),lb,ub,lwd=2)]
  # Legend.
  DT[,legend('topright',legend=c(sprintf("Overall mean %1.2f",mean(get(yvar),na.rm=T)),sprintf("Overall median %1.2f",median(get(yvar),na.rm=T))),cex=1.5,pch=NULL,bty='n')]
  return(bias)
}


####################### FIGURE 2 #####################  
# Scatter plot of subjective prob on observed rate.
pdf(file=file.path(.fig.path,sprintf("SubjectiveProbOnVoteSplitGivenChars%s.pdf",postpend)),width=9,height=6)
par.old = par(mar=c(5.1,4.1,1.1,1.1))

binned = DT[,list(mn.trmp=mean(posterior.trump,na.rm=T),mdn.trmp=median(posterior.trump,na.rm=T),.N),keyby=round(100*trump.2p,0)]
# Plot.
binned[,plot(x=round,y=mn.trmp,ann=F,axes=F,type='n',ylim=range(c(mn.trmp,mdn.trmp)))]
grid()
axis(2,las=2);axis(1)
title(xlab="Observed ANES Trump rate given characteristics",ylab="Subject probability Trump given characteristics")
# Points through means
binned[,points(x=round,y=mn.trmp,pch=19,cex=1.5)]
# Loess through means
binned[,lines(loess.smooth(x=round,y=mn.trmp),lwd=4,col=3)]

# Regression slope:
legend("bottomright",legend=sprintf("Regression slope %1.2f",binned[,coef(lm(mn.trmp~round))[2]]),pch=NA,bg='white')
rm(binned)
par(par.old)
dev.off()

# Regression slope for relationship.
print(summary(DT[,lm(posterior.trump ~ I(100*trump.2p))]))


################# TABLE A1 (NEEDED FOR FIGURE 3) ###################  
library(rockchalk)
varLabels = list("treat.gender"="Gender","treat.mip"="Most important problem","treat.pid"="Party identification","treat.income"="Income","treat.race"="Race","treat.state"="State",prior="Prior","l.prior"="Logit prior")
out <- outreg(list(
"Basic"=DT[,lm(posterior ~ -1 + prior + treat.pid + treat.mip + treat.race + treat.gender + treat.income + treat.state)],
"Prior less than 33"=DT[prior < 33,lm(posterior ~ -1 + prior + treat.pid + treat.mip + treat.race + treat.gender + treat.income + treat.state)],
"Prior 33 to 67"=DT[prior >= 33 & prior < 67,lm(posterior ~ -1 + prior + treat.pid + treat.mip + treat.race + treat.gender + treat.income + treat.state)],
"Prior greater than 67"=DT[prior >= 67,lm(posterior ~ -1 + prior + treat.pid + treat.mip + treat.race + treat.gender + treat.income + treat.state)],
# By payment condition.
"Paid per round"=DT[bonus=="$0.10",lm(posterior ~ -1 + prior + treat.pid + treat.mip + treat.race + treat.gender + treat.income + treat.state)],
"Flat payment"=DT[bonus=="one point",lm(posterior ~ -1 + prior + treat.pid + treat.mip + treat.race + treat.gender + treat.income + treat.state)]), 
type="latex", varLabels=varLabels, centering="dcolumn", print.results=F,digits=2,alpha=c(0.05))
# Remove the "Estimate" and "S.E." rows.
out <- out[regexpr("Estimate",out) == -1]; out <- out[regexpr("S\\.E\\.",out) == -1]
# Write to file here.
write(out,file=file.path(.tab.path,sprintf("LearningTruth%s.tex",postpend)))


####################### FIGURE 3 #####################  
# Coefficient plot of estimated informativeness of each input value.
pdf(file.path(.fig.path,sprintf("InformativenessByInfoValue%s.pdf",postpend)),width=8,height=6)
par.old = par(mar=c(3.1,6.65,0.6,1.1))
# Coefficients.
coefs = .mkCoefs(DT[,lm(posterior ~ -1 + prior + full.treat)],regex="full.treat.*: |full.treat")
# Break in two frames.
.coefPlot(coefs[seq(floor(.N/2)+1,.N),],ylab="",xlab="Informativeness towards truth",ynames="yname",cex.axis.y=.9)
.coefPlot(coefs[seq(1,floor(.N/2)),],ylab="",xlab="Informativeness towards truth",ynames="yname",cex.axis.y=.9)
par(par.old)
dev.off()


########################## FIGURES 4 and A9 #################
# Accuracy by number of shared characteristics.

# Cumulative number of shared characteristics.
setkey(DT,our.id,round,elicitation.num)
DT[is.na(in.group), in.group := 0]
DT[!is.na(in.group), cumsum.in.group := cumsum(in.group), by=c("our.id","round")]
DT[,fac.cumsum.in.group := as.character(cumsum.in.group)]
DT[fac.cumsum.in.group %in% c("3","4"),fac.cumsum.in.group := "3-4"]
DT[,fac.cumsum.in.group := factor(fac.cumsum.in.group)]
# Aggregate mean accuracy to elicitation number and
# number shared chars.
accsc = DT[!is.na(prob.in.ci),.(in.ci=mean(prob.in.ci),num.elicits=.N),keyby=c("elicitation.num","fac.cumsum.in.group")]

# Append aggregated accuracy across elicitation.nums.
accsc = rbind(accsc,DT[!is.na(prob.in.ci),.(in.ci=mean(prob.in.ci),num.elicits=.N),keyby=c("fac.cumsum.in.group")],fill=T)
accsc[,elicitation.num := as.character(elicitation.num)]
accsc[is.na(elicitation.num), elicitation.num := "Overall"]
# Calculate error bar bounds.
accsc[,ub := .CI95(in.ci,num.elicits)]
accsc[,lb := .CI95(in.ci,num.elicits,upper=F)]

# Add in empty cells for shared chars in early
# elicitations. For plot purposes.
joe = merge(accsc[,unique(elicitation.num)], accsc[,unique(fac.cumsum.in.group)])
setnames(joe,c("elicitation.num","fac.cumsum.in.group"))
accsc = merge(accsc,joe,by=c("elicitation.num","fac.cumsum.in.group"),all=T)
rm(joe)

# Average distance between estimate and target
# by number of shared characteristics and elicitation
# number.
DT[elicitation.num > 1,elicit.cumsum := factor(sprintf("%s-%s",elicitation.num, gsub("3-4","3/4",fac.cumsum.in.group)), levels= c("2-0", "2-1", "3-0", "3-1", "3-2", "4-0", "4-1", "4-2", "4-3/4", "5-0", "5-1", "5-2", "5-3/4"))]
pdf(file=file.path(.fig.path,sprintf("BiasBySharedChars%s.pdf",postpend)),width=10,height=7)
par.old = par(mar=c(4.1,4.1,0.6,0.6))
xvar = "elicit.cumsum"

# Any estimate.
b = .biasPlot(DT[!is.na(get(xvar)),],yvar="dist.from.share",xvar=xvar,xlab="Elicitation number-Number shared characteristics",add.lines=F)
# Add lines connecting within elicitation.num.
b[,e.num := substr(elicit.cumsum,1,1)]
b[input != "Overall",lines(x=as.numeric(get(xvar)),y=dist.mean,col='gray'),by=c("input","e.num")]
b[input == "Overall",lines(x=as.numeric(get(xvar)),y=dist.mean,lwd=3),by=e.num]
# Re-do points on top.
b[input == "Overall",points(x=as.numeric(get(xvar)),y=dist.mean,pch=19,cex=1.5)]

# Only estimates outside of CI.
b = .biasPlot(DT[!is.na(get(xvar)),],yvar="dist.from.share.notin.ci",xvar=xvar,xlab="Elicitation number-Number shared characteristics",add.lines=F)
# Add lines connecting within elicitation.num.
b[,e.num := substr(elicit.cumsum,1,1)]
b[input != "Overall",lines(x=as.numeric(get(xvar)),y=dist.mean,col='gray'),by=c("input","e.num")]
b[input == "Overall",lines(x=as.numeric(get(xvar)),y=dist.mean,lwd=3),by=e.num]

#
# Average distance between estimate and target
# by number of shared characteristics and subject vote.
# 
DT[elicitation.num > 1,vote.16.cumsum := factor(sprintf("%s\n%s",vote.16, gsub("3-4","3/4",fac.cumsum.in.group)), levels= c("Clinton\n0", "Clinton\n1", "Clinton\n2", "Clinton\n3/4","Trump\n0", "Trump\n1", "Trump\n2", "Trump\n3/4","Other\n0", "Other\n1", "Other\n2", NA))]
DT[,table(vote.16.cumsum)]
xvar = "vote.16.cumsum"

# Any estimate.
b = .biasPlot(DT[!is.na(get(xvar)),],yvar="dist.from.share",xvar=xvar,xlab="Subject vote and Number shared characteristics",add.lines=F,tick=F)
# Add lines connecting within subject vote.
b[,vote.16 := ifelse(grepl("Clinton",vote.16.cumsum),"Clinton", ifelse(grepl("Trump",vote.16.cumsum),"Trump", "Other"))]
b[input != "Overall",lines(x=as.numeric(get(xvar)),y=dist.mean,col='gray'),by=c("input","vote.16")]
# Re-do points on top.
b[input == "Overall",lines(x=as.numeric(get(xvar)),y=dist.mean,lwd=3),by=vote.16]
b[input == "Overall",points(x=as.numeric(get(xvar)),y=dist.mean,pch=19,cex=1.5)]

par(par.old)
dev.off()

##################### FIGURE A1 ####################  
# Bar chart of percent subjective probabilities w/in
# 95\% confidence interval of vote share given other's
# characteristics, by elicitation number.
  
# Aggregate to elicitation number.
in.ci = DT[!is.na(prob.in.ci),.(in.ci=mean(prob.in.ci),num.elicits=.N),keyby=c("round","elicitation.num")]
# Append overall.
in.ci = rbind(in.ci,DT[!is.na(prob.in.ci),.(in.ci=mean(prob.in.ci),num.elicits=.N),by=c("elicitation.num")],fill=T)
in.ci[is.na(round), round := 5]
# Calculate error bar bounds.
in.ci[,ub := .CI95(in.ci,num.elicits)]
in.ci[,lb := .CI95(in.ci,num.elicits,upper=F)]

# Bar plot, grouped by round and then elicitation number.
pdf(file=file.path(.fig.path,sprintf("PctInCIByRound%s.pdf",postpend)),width=10,height=6)
par.old = par(mar=c(4.1,4.1,1.1,1.1))
z = in.ci[,barplot(in.ci~elicitation.num+round,beside=T,xlab="Round and elicitation number",ylab="Proportion within confidence interval",axes=F,names.arg=c("1","2","3","4","Overall"),ylim=range(round(c(ub,lb,0),3)))]
axis(2,las=2)
# Bar plot error bars.
in.ci[,.errorBars(z=z,ub,lb)]
par(par.old)
dev.off()

  
############## FIGURES A2 and A6 ############
# Bar chart of average movement towards truth by
# prior and characteristic.
# Learning by prior, elicitation, characteristic.
DT[,prior.cat := cut(prior,breaks=seq(0,100,10),include.lower=T)] # tenths.
movement2 = DT[!is.na(prior) & !is.na(posterior),list(dist.toward.true.mn=mean(dist.toward.true),dist.toward.true.sd=sd(dist.toward.true),.N),keyby=c("input","prior.cat","elicitation.num")]
# Create confidence interval.
movement2[,mn.se := dist.toward.true.sd/sqrt(N)]
movement2[,lb := dist.toward.true.mn - 1.96*mn.se]
movement2[,ub := dist.toward.true.mn + 1.96*mn.se]
# Overall means.
move3 = DT[!is.na(prior) & !is.na(posterior),list(dist.toward.true.mn=mean(dist.toward.true),dist.toward.true.sd=sd(dist.toward.true),.N),keyby=c("input")]
move3[,rev.sort := -1*dist.toward.true.mn]
setkey(move3,rev.sort)
# Create a factor variable input sorted descending
# on overall mean toward T.
movement2[,input.fac := factor(toupper(input),levels=toupper(move3[,input]))]
move3[,input.fac := factor(toupper(input),levels=toupper(move3[,input]))]

# Learning by prior and characteristic only.
movement3 = DT[!is.na(prior) & !is.na(posterior),list(dist.toward.true.mn=mean(dist.toward.true),dist.toward.true.sd=sd(dist.toward.true),.N), keyby=c("input","prior.cat")]
# Create confidence interval.
movement3[,mn.se := dist.toward.true.sd/sqrt(N)]
movement3[,lb := dist.toward.true.mn - 1.96*mn.se]
movement3[,ub := dist.toward.true.mn + 1.96*mn.se]
# Create a factor variable input sorted descending
# on overall mean toward T.
movement3[,input.fac := factor(toupper(input),levels=toupper(move3[,input]))]

pdf(file=file.path(.fig.path,sprintf("MeanMoveTowardsTruth%s.pdf",postpend)),width=10,height=7)
par.old = par(mar=c(2.1,4.1,0.3,0.1))
.yvar = "dist.toward.true.mn"
# One version without error bars.
setkey(movement3,prior.cat,input.fac) # sorted on p.true
z = movement3[, barplot(get(.yvar)~input.fac+prior.cat,beside=T,xlab="Prior belief and elicitation number",ylab="",axes=F,legend.text=F)]
axis(2,las=2)

# Second version with error bars.
z = movement3[, barplot(get(.yvar)~input.fac+prior.cat,beside=T,xlab="Prior belief and elicitation number",ylab="",axes=F,legend.text=F)] 
axis(2,las=2)
# Bar plot error bars.
movement3[,.errorBars(z=z,ub,lb,lwd=1)]

par(par.old)
dev.off()


####################### FIGURES A3 and A8 ###############  
# Line plot of average distance between estimate 
# and vote share (y-axis) given prior beliefs (binned on
# x-axis). Then, separate lines for each characteristic
# for rounds having just observed that characteristic.
DT[,prior.cat := NULL]
# Prior into 20ths.
DT[!is.na(trump.2p),prior.cat := cut(prior,breaks=seq(0,100,5),include.lower=T)] # twentieths.
  
pdf(file=file.path(.fig.path,sprintf("BiasByPrior%s.pdf",postpend)),width=10,height=7)
par.old = par(mar=c(4.1,4.1,0.6,0.6))
# Average distance.
# Any estimate.
b = .biasPlot(DT[!is.na(prior.cat),],"dist.from.share")
# Only estimates outside of CI.
b = .biasPlot(DT[!is.na(prior.cat),], "dist.from.share.notin.ci")

# Average absolute distance.
# Problem: This includes all the 
#b = .biasPlot(DT[!is.na(prior.cat),], "abs.dist.from.share", ylab="Absolute value average estimate minus ANES share")
par(par.old)
dev.off()


#
# Table of bias by subject 2016 vote. Referenced in text.
#
cat("\n\nTabulation of bias by subject 2016 vote:\n")
.biasTab <- function(DT,yvar="dist.from.share",xvar="vote.16") {
  # Create the bias table.
  b = DT[, list("dist.mean"=mean(get(yvar),na.rm=T),"dist.median"=median(get(yvar),na.rm=T), "sd"=sd(get(yvar),na.rm=T),.N), keyby=xvar]
  b[,SE := sd/sqrt(N)]
  return(b)
}
print(.biasTab(DT))

######################## FIGURE A4 ##############
# Distribution of changes in beliefs by input type.
pdf(file.path(.fig.path,sprintf("InitialBeliefs%s.pdf",postpend)),width=10,height=6)
par.old = par(ask=F,mfrow=c(1,1),oma=c(0,0,0,0),mar=c(4.1,4.1,1.1,1.1),cex.axis=.95)
# Breaks.
the.breaks = c(seq(0,40,10),45,49,51,55,seq(60,100,10))
# Break into categories.
.tmp = cut(DT[elicitation.num == 1,posterior.trump],breaks=the.breaks)
DT[elicitation.num == 1,post.trump.cats := .tmp]
# Overall density.
tab = DT[elicitation.num == 1,prop.table(table(post.trump.cats))]
plot(tab,type="l",xlab="Initial probability Trump",ylab="Proportion",lwd=3,ylim=c(0,0.55))
# Densities by subject vote choice.
gcolors = rev(gray(seq(.45,.65,length.out=3)))
for (trt in (.inputs = c("Clinton","Trump","Other"))) {
  .tab = DT[vote.16 == trt & elicitation.num == 1,prop.table(table(post.trump.cats))]
  lines(.tab,type="l",lwd=3,col=gcolors[which(trt == .inputs)])
}
legend("topright",lwd=3,col=c("black",gcolors),legend=c("OVERALL",toupper(.inputs)))
# Mover ovarll on top again.
lines(tab,type="l",lwd=3)

# Separate plot initial beliefs by round.
plot(tab,type="l",xlab="Initial probability Trump",ylab="Proportion",lwd=3,ylim=c(0,0.55))
gcolors = rev(gray(seq(.45,.65,length.out=4)))
for (rnd in (.inputs = c("1","2","3","4"))) {
  .tab = DT[round == rnd & elicitation.num == 1,prop.table(table(post.trump.cats))]
  lines(.tab,type="l",lwd=3,col=gcolors[which(rnd == .inputs)])
}
legend("topright",lwd=3,col=c("black",gcolors),legend=c("OVERALL",sprintf("ROUND %s",1:4)))
# Move overall on top again.
lines(tab,type="l",lwd=3)

par(par.old)
dev.off()

################### FIGURE A5 ###############  
# Histogram of ANES vote splits for characteristic 
# combinations delivered to subjects.
# Page 2 is among round 2 and later probs elicted from respondents.
pdf(file=file.path(.fig.path,sprintf("Histogram-VoteSplitsGivenChars%s.pdf",postpend)),width=7,height=6)
par.old = par(mfrow=c(2,1),mar=c(4.1,4.1,2.1,0.6))
# ANES.
trump.shares = fread("TrumpShareGivenChars.csv")
x = trump.shares[,hist(100*trump.2p,breaks=30,col='gray',xlab="",main="Distribution of vote splits\namong characteristic sets seen by subjects")]
title(xlab="Trump two-party share",line=2)
# Study subjects. Only for elicitations after first.
DT[elicitation.num > 1,hist(posterior.trump,breaks=x$breaks,col='gray',xlab="",main="Distribution of probabilities returned by subjects")]
title(xlab="Elicited probability Trump",line=2)
par(par.old)
dev.off()
  
# Trump 2p share for females from Texas for paper text.
cat("\n++++++++++++++++\nTrump 2-party share for females from Texas, ANES:\n")
print(trump.shares[the.values == "Female|::|Texas",])


################### FIGURE A7 ################
# Bar chart of average movement towards truth by
# prior, elicitation number, and characteristic.
pdf(file=file.path(.fig.path,sprintf("MeanMoveTowardsTruth-ByElicitation%s.pdf",postpend)),width=10,height=2)
par.old = par(mar=c(2.1,4.1,0.3,0.1))
.yvar = "dist.toward.true.mn"
# Loop over elicitation.num and make bar chart.
# One version without error bars.
setkey(movement2,prior.cat,input.fac) # sorted on p.true
for (i in 2:5) { 
  z = movement2[elicitation.num == i, barplot(get(.yvar)~input.fac+prior.cat,beside=T,xlab="Prior belief and elicitation number",ylab="",axes=F,legend.text=F)]
  axis(2,las=2)
}
# Second version with error bars.
for (i in 2:5) { 
  z = movement2[elicitation.num == i, barplot(get(.yvar)~input.fac+prior.cat,beside=T,xlab="Prior belief and elicitation number",ylab="",axes=F,legend.text=F)]
  
  axis(2,las=2)
  # Bar plot error bars.
  movement2[elicitation.num == i,.errorBars(z=z,ub,lb,lwd=1)]
}
par(par.old)
dev.off()




# Sample size and item non-response.
n.respondent = DT[,length(unique(our.id))]
cat("\n\nNumber of respondents in sample:",n.respondent,"\n")
cat("Number of prior beliefs:",DT[!is.na(prior),.N],"\n")
cat("Expected number of prior beliefs:",4*4*n.respondent,"\n")
cat("Number of posterior beliefs:",DT[!is.na(posterior),.N],"\n")
cat("Expected number of prior beliefs:",5*4*n.respondent,"\n\n")

